function [f,f_n,Aineq_d,Aineq_i,bineq,lb,ub,ctype,dr_hh,dr,Xu,gu_p,gr3,gr6,gr7,Xr3,Xr6,Xr7,Yr3,Yr6,Yr7,Dr3,Dr6,Dr7] = generate_cplex_inputs_weight(bsperm,g,X,bsi,k,gn,Xn,g_bywave,X_bywave,Y_bywave,D_bywave,opower_paper_wave,wave,sample_wave,target_wave)
% This function generates inputs required for the CPLEX solver, after
% applying density ratio

% Input: 
% (1) g: treatment effect calculated using IPW 
% (2) X: set of covariates 
% (3) k: number of covariates 
% (4) g_byset, x_byset, Y_byset, D_byset, opower_byset, set_byset: key 
%                   inputs by wave
% (5) opower_paper_wave: opower_wave in the original sample order
% (6) wave: wave=3,6,or7, indicates need to reweight to a particular wave, 
% not the pooled wave
% (7) sample_wave: historical wave 
% (8) target_wave: future wave 

% Output list is too long, includes all wave-specific cplex inputs.

%% Generate permutated sample (for CI calculation)
    g_perm = g(bsperm,:);
    X_perm = X(bsperm,:);
    opower_paper_wave_perm = opower_paper_wave(bsperm);
    % Subtract Wn from Wbs
    if (bsi>0)
        g_perm = [g_perm; -gn]; 
        X_perm = [X_perm; Xn];
    end
%% Generate gu and nXu for unique households (same X) 
    [Xr, Ind] = sortrows(X_perm);
    gr = g_perm(Ind);
    Xu = unique(Xr,'rows','stable');
    nu  = size(Xu,1); % number of unique covariate vectors
    gu = zeros(nu,1);
    
    Xr3=Xr(opower_paper_wave_perm(Ind)==3,:);
    Xr6=Xr(opower_paper_wave_perm(Ind)==6,:);
    Xr7=Xr(opower_paper_wave_perm(Ind)==7,:);
    n=length(g); n3=size(Xr3,1); n6=size(Xr6,1); n7=size(Xr7,1);
    % stack three waves together for opower_paper_wave, X_bywave,
    % g_bywave, then sort within wave 
    opower_bywave=[repelem(3,size(Xr3,1))';repelem(6,size(Xr6,1))';repelem(7,size(Xr7,1))'];
    [Xr_bywave,Ind_bywave]=sortrows(X_bywave);
    gr_bywave=g_bywave(Ind_bywave);
    Yr_bywave=Y_bywave(Ind_bywave);
    Dr_bywave=D_bywave(Ind_bywave);

    gr3=gr_bywave(opower_bywave(Ind_bywave)==3,:);
    gr6=gr_bywave(opower_bywave(Ind_bywave)==6,:);
    gr7=gr_bywave(opower_bywave(Ind_bywave)==7,:);
    
    Yr3=Yr_bywave(opower_bywave(Ind_bywave)==3,:);
    Yr6=Yr_bywave(opower_bywave(Ind_bywave)==6,:);
    Yr7=Yr_bywave(opower_bywave(Ind_bywave)==7,:);
    
    Dr3=Dr_bywave(opower_bywave(Ind_bywave)==3,:);
    Dr6=Dr_bywave(opower_bywave(Ind_bywave)==6,:);
    Dr7=Dr_bywave(opower_bywave(Ind_bywave)==7,:);
    
    if wave==0 % if in the pooled sample, all waves adjust to the overall X distribution
        [gu3,nXu3] = empirical_density_generate_gu_cubic(Xr3,Xu,gr3);
        [gu6,nXu6] = empirical_density_generate_gu_cubic(Xr6,Xu,gr6);
        [gu7,nXu7] = empirical_density_generate_gu_cubic(Xr7,Xu,gr7);
        [gu_original,nXu] = empirical_density_generate_gu_cubic(Xr,Xu,gr);
        
        dr3=((nXu3+nXu6+nXu7)./n)./(nXu3./n3); 
        dr6=((nXu3+nXu6+nXu7)./n)./(nXu6./n6); 
        dr7=((nXu3+nXu6+nXu7)./n)./(nXu7./n7);
        dr3(isinf(dr3))=0;
        dr6(isinf(dr6))=0;
        dr7(isinf(dr7))=0;
        
        gu3_p=gu3.*dr3;
        gu6_p=gu6.*dr6;
        gu7_p=gu7.*dr7;
        % The pooled gu is the sum of wave-specific gu. 
        gu_p=gu3_p+gu6_p+gu7_p;
   
    else % if wave analysis, then weight to a particular wave
     if sample_wave==3 & target_wave==6 
        [gu6,nXu6,ic6,wave6_in_pool] = empirical_density_generate_gu_cubic(Xr6,Xu,gr6);
        [gu3,nXu3,ic3,wave3_in_pool] = empirical_density_generate_gu_cubic(Xr3,Xu,gr3);
        dr=((nXu6)./n6)./(nXu3./n3); 
        dr(isinf(dr))=0;
        dr(isnan(dr))=0;
        %dr_sample_wave_target_wave
        gu_p=gu3.*dr;
        % Re-weight gu of sample population to get the gu for the target
        % population, then calculate the ewm rules 
        
      elseif sample_wave==3 & target_wave==7 
        [gu7,nXu7,ic7,wave7_in_pool] = empirical_density_generate_gu_cubic(Xr7,Xu,gr7);
        [gu3,nXu3,ic3,wave3_in_pool] = empirical_density_generate_gu_cubic(Xr3,Xu,gr3);
        dr=((nXu7)./n7)./(nXu3./n3); 
        dr(isinf(dr))=0;
        dr(isnan(dr))=0;
        gu_p=gu3.*dr;
        
      elseif sample_wave==6 & target_wave==3  
        [gu3,nXu3,ic3,wave3_in_pool] = empirical_density_generate_gu_cubic(Xr3,Xu,gr3);
        [gu6,nXu6,ic6,wave6_in_pool] = empirical_density_generate_gu_cubic(Xr6,Xu,gr6);
        dr=((nXu3)./n3)./(nXu6./n6); 
        dr(isinf(dr))=0;
        dr(isnan(dr))=0;
        gu_p=gu6.*dr;         
      elseif sample_wave==6 & target_wave==7 
        [gu7,nXu7,ic7,wave7_in_pool] = empirical_density_generate_gu_cubic(Xr7,Xu,gr7);
        [gu6,nXu6,ic6,wave6_in_pool] = empirical_density_generate_gu_cubic(Xr6,Xu,gr6);
        dr=((nXu7)./n7)./(nXu6./n6); 
        dr(isinf(dr))=0;
        dr(isnan(dr))=0;
        gu_p=gu6.*dr;
      elseif sample_wave==7 & target_wave==3 
        [gu3,nXu3,ic3,wave3_in_pool] = empirical_density_generate_gu_cubic(Xr3,Xu,gr3);
        [gu7,nXu7,ic7,wave7_in_pool] = empirical_density_generate_gu_cubic(Xr7,Xu,gr7);
        dr=((nXu3)./n3)./(nXu7./n7); 
        dr(isinf(dr))=0;
        dr(isnan(dr))=0;
        gu_p=gu7.*dr;
      elseif sample_wave==7 & target_wave==6 
        [gu6,nXu6,ic6,wave6_in_pool] = empirical_density_generate_gu_cubic(Xr6,Xu,gr6);
        [gu7,nXu7,ic7,wave7_in_pool] = empirical_density_generate_gu_cubic(Xr7,Xu,gr7);
        dr=((nXu7)./n7)./(nXu6./n6);
        dr(isinf(dr))=0;
        dr(isnan(dr))=0;
        gu_p=gu7.*dr;
      end  
    end
%% use gu and Xu to set up the cplex inputs matrix
    % Add explicit monotonicity constraints
    % Xu is ordered by increasing covariate (Xu(:,2))
    % and then by increasing usage (Xu(:,5))
    % For each level of covariate we impose treatment set inclusion
    sameinc = (Xu(1:end-1,2)==Xu(2:end,2));
    % decreasing in pre-treatment consumption
    Mineq_d = [diag(sameinc) zeros(nu-1,1)] + [zeros(nu-1,1) diag(-sameinc)];
    % increasing in pre-treatment consumption
    Mineq_i = [diag(-sameinc) zeros(nu-1,1)] + [zeros(nu-1,1) diag(sameinc)];

    f = [zeros(k,1); gu_p]; % objective function coefficients
    f_n = [zeros(k,1); -gu_p]; % Reverse objective function coefficients
    B = 1; % bounds on coefficients
    C = B*sum(abs(Xu),2); % maximum values of x'beta
    minmargin = max(1,C)*(1e-8); 
    Aineq_d = [[-Xu diag(C)]; [Xu -diag(C)]; [zeros(nu-1,k) Mineq_d]]; 
    Aineq_i = [[-Xu diag(C)]; [Xu -diag(C)]; [zeros(nu-1,k) Mineq_i]];
    bineq = [[C-minmargin];[-minmargin]; minmargin(1:nu-1,:)];
    lb = [-B*ones(k,1); zeros(nu,1)];
    ub = [ B*ones(k,1);  ones(nu,1)];
    
    % Variable type string
    ctype = strcat(repmat('C',1,k),repmat('B',1,nu));
    
%% Get household level dr (as opposed to unique household level dr used to generate the cplex input matrix above)
switch sample_wave
    case 3
        dr3 = dr(wave3_in_pool);
        dr_hh=dr3(ic3);
    case 6
        dr6 = dr(wave6_in_pool);
        dr_hh=dr6(ic6);
    case 7
        dr7 = dr(wave7_in_pool);
        dr_hh=dr7(ic7);
end

end 

